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A two-offspring branching annihilating random walk model, with finite reaction rates, is studied 
in one-dimension. The model exhibits a transition from an active to an absorbing phase, expected 
to belong to the DP2 universality class embracing systems that possess two symmetric absorbing 
states, which in one-dimensional systems, is in many cases equivalent to parity conservation. The 
' f-H ' ' phase transition is studied analytically through a mean-field like modification of the so-called parity 

, interval method. The original method of parity intervals allows for an exact analysis of the diffusion- 

■ controlled limit of infinite reaction rate, where there is no active phase and hence no phase transition. 

For finite rates, we obtain a surprisingly good description of the transition which compares favorably 
with the outcome of Monte Carlo simulations. This provides one of the first analytical attempts to 
deal with the broadly studied DP2 universality class. 
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I. INTRODUCTION 



S 

Phase transitions occurring away from thermodynamical equilibrium constitute one of the most challenging topics 
in statistical physics. They appear in a host of physical systems as well as in many models in biology, chemistry, 
O , sociology, etc. Given the lack of a general theory of non-equilibrium systems, we are still living a taxonomic era in this 
field: it would be highly desirable to reach a complete classification of the known phase transitions into universality 
| classes, as a preliminary step to their full categorization, and identification of relevant features. 
^ ■ After more than twenty years of study, it has become clear that the degree of universality is much narrower away 
from equilibrium than it is in equilibrium transitions. Certainly, general properties such as symmetries, conservation 
laws, dimensionalities, etc., play a key role, as they do in equilibrium. But some other ingredients, such as microscopic 
dynamical details, hard core interactions, and the type of updating may, in some cases, influence the emerging critical 
behavior of non-equilibrium systems, making the task of theoreticians both stimulating and difficult. 
• One of the most robust and best studied non-equilibrium classes of transitions is directed percolation (DP). It 
includes an amazing variety of models and systems exhibiting a transition from an active to an absorbing phase 
"p-; ■ [1, 2, 3, 4, 5, 6], and is well characterized at a field-theoretical level by Reggeon field theory (RFT) [7, 8]. 

DP is so robust, that identifying and classifying the nature of perturbations able to drive absorbing phase transitions 
away from this class has become a challenging task. Probably the best known instance of this is the DPI class, where 
the presence of two perfectly symmetric absorbing states (Z2 symmetry) is the main feature responsible for non-DP 
scaling [9]. This class includes a cellular automaton introduced by Grassberger et a!., [10], interacting monomer-dimer 
^ models [11], non-equilibrium Ising models [12, 13], monomer-monomer surface reaction models [14], ^-symmetric 
generalizations of DP [15] and of the contact process [16], branching annihilating random walks with conserved parity 
[17, 18, 19], and it has also been related to generalized versions of the Voter model [20, 21]. Extensive numerical 
simulations led to conjecture rational exponent values for the DP2 class [22], though this was later disproved by more 
exhaustive simulations [4]. 

Although the existence and robustness of the DP2 class arc well established from a numerical view-point, a solid 
theoretical understanding is still missing. Bold attempts to write down and renormalize a field theory suitable for the 
DP2 class have been performed, but the results are not as satisfactory as they are for RFT [23, 24]. In particular, the 
renormalization is based on a clever, but somehow uncontrolled expansion around two-different critical dimensions. 
As we shall illustrate below, straightforward mean-field approaches (cluster approximations) fail to reproduce a phase 
transition if too small clusters are considered, and one has to resort to large-clusters, which make the calculation 
complicated and not very accurate, i.e. the convergence of the series is very slow (although results can be improved 
if combined with coherent- anomaly methods [13]). 

It is the purpose of this paper to shed some light on these issues, by examining models of Branching Annihilating 
Walks with two offspring (2-BAW) in one dimension: the 2-BAW with finite reaction rate exhibits a phase transition 
in the DP2 class [25]. Our approach is based on a mean-field modification of the method of parity intervals, originally 
introduced for the study of annihilation reactions, A + A ^ [26, 27, 28, 29]. A similar approach has been previously 
employed in conjunction with the method of empty intervals [26, 30, 31] for the analysis of other intractable models 
that do not conserve parity [32]. See [27, 28] for a more detailed introduction to this method and its applications to 
different models. 
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The paper is organized as follows. The 2-BAW model is described in Section II, where we also perform simple 
(up to two sites) cluster approximations (which fail to capture the transition). In Section III we present the exact 
solution of the 2-BAW model in the limit of infinite reaction rate, using the method of parity-intervals. Although 
there is no phase transition in this limit, the analysis serves as a basis for the mean- field like approximation presented 
in Section IV, for the relevant case of finite reaction rates. There we derive the approximate steady-state solution for 
generic parameter values, and compare it with the outcome of Monte Carlo computer simulations. We conclude with 
a critical discussion of our results and further developments, in Section V. 

II. THE MODEL AND CLUSTER APPROXIMATIONS 

The model is defined as follows. Each site of a one-dimensional lattice is either empty or singly occupied. The lattice 
is updated asynchronously: an occupied site is randomly chosen, and it is tried for diffusion, at rate T (probability 
r/(r + 0)), or branching, at rate £1 (probability 0/(r + f2)); time is increased by l/N, where N is the number of 
occupied (active) sites. In a diffusion attempt the particle is moved to one of its two nearest neighbors (target sites), 
with equal probabilities. If the target site is occupied, annihilation results with probability r: the move is effected and 
both particles, the diffuser and the target, are removed from the lattice. The move is rejected with probability 1 — r, 
and the lattice state remains unchanged. In a branching attempt the particle gives birth to two new particles onto 
the nearest neighbor sites (target sites). If either, or both of the target sites is occupied, annihilation takes place with 
probability r: branching is effected, and the occupied target site(s) become empty. Once again, the move is rejected 
with probability 1 — r, and the lattice state remains unchanged. 

In the diffusion-controlled limit of infinite reaction rate (r = 1), the 2-BAW model is known to possess only a steady 
absorbing phase [33]. However, for finite rates (r < 1), a transition belonging to the DP2 class, from an absorbing 
state, at r < r c , to an active phase, at r > r c , is found in numerical simulations [25]. We now demonstrate that simple 
mean-field approximations fail to capture this transition. 

The simplest conceivable mean-field theory is that of the one-site approximation, obtained by neglecting all cor- 
relations between the states of different sites. For example, if the probability of one site being occupied is equal 
to the concentration c, then the probability of finding one occupied site followed immediately by an empty site is 
c(l — c). There are only three events that lead to a change in particle concentration: (i) annihilation of two particles 
by diffusion, •• — > oo; (ii) creation of two particles via birth, o»o — > and (iii) annihilation of two particles due to 
birth onto previously occupied sites, • • • — > o • o. Note that the birth process o o (and its mirror-symmetric 

image) does not alter the concentration. However, processes (i) and (iii) occur with restricted probability r. Thus 

-^c = -2rc 2 + 2c(l - cf - 2rc 3 , (1) 

where we took F = O = 1. The first, second and third terms on the r.h.s. correspond to (i), (ii) and (iii), respectively. 

For the steady state, we set dc/dt = and choose the stable root (Fig. 1), c s = (2 + r — \/8r + r 2 )/(2 — 2r). We 
conclude that according to the one-site approximation the system always evolves to an active phase, with 1/3 < c s < 1 
(Fig.2). 

The two-site approximation provides the next level of complexity. Let p be the conditional probability for a site to 
be occupied, given that the adjacent site (to its left) is occupied. Let q be the conditional probability for a site to 
be occupied, given that the site to its left is empty. Thus, the probability of the events ••, »o, o», oo is cp, c(l — p), 
(1 — c)q, (1 — c)(l — q), respectively. Note that since »o and o» are equally likely, we have, 

c(l-p) = (l-c)q. (2) 

(The same relation may be derived from the fact that Pr(») = Pr(»») + Pr(o»).) Eq. (1) is now rewritten as 

-|c = -2rcp + 2(1 - c)q{\ - p) - 2rcp 2 , (3) 

where the terms on the r.h.s. correspond to the same events as before. Setting dc/dt = and using the relation (2), 
we obtain, for the steady state, 

= -2rcp + 2c(l - pf - 2rcp 2 . 

We conclude that either c s = or p s = (2 + r — V8r + r 2 )/(2 — 2r). In the latter case, an additional evolution 
equation, for the event ••, provides the missing value of c s : c s > for all r. The result agrees closely with that of the 
one-site approximation, with deviations not larger than 20%. Thus, also the two-site approximation fails to predict 
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FIG. 1: One-site cluster approximation. Flow diagram, according to Eq. (1), with dc/dt — 0. This shows that the root c s = 
is unstable, while the root c s > is stable. Results from the two-site approximation are very similar. The stationary solution 
as a function of r is plotted in Fig. 2 



the transition. Considering larger clusters, a phase transition can be generated [13, 34], but a large set of coupled 
equations (which enlarges with cluster-size) has to be solved, the accuracy is not good, and the results are expected 
to be valid only up to the cluster-size. In the next section we present an alternative method intended to overcome 
these difficulties. 



III. PARITY INTERVALS AND EXACT SOLUTION FOR r = 1 



We now turn to a different approach, that of the method of parity intervals [26, 27, 28, 29]. We first present the 
exactly soluble case of r = 1, for which there is no transition. The exact approach followed here serves as a basis for 
a mean-field like approximation, for the more interesting case of r < 1 — an approximation that does capture the 
transition and reproduces kinetic details surprisingly well. 

Let G n (t) be the probability that (in an homogeneous system) an arbitrary segment of n consecutive sites contains 
an even number of particles at time t. A site can be either empty or occupied by a single particle, so the probability 
that a site is occupied, i.e., the particle density, is 

p(t) = l-G 1 (t). (4) 

Since the dynamic rules of the two-offspring BAW conserve parity, the only way that G n might change is when: 
(i) particles at the edge of the segment hop outside or branch, or (ii) particles just outside of the segment hop inside 
or branch. Let i i and ^^^h represent segments of even and odd number of particles, respectively. Then, the 

changes of G n can be described schematically by the events: 




n — 1 ' n * 

] • — I 



(5) 



Arrows in the first four terms indicate the hopping of a particle to the left or right. In the last four terms, the arrows 
indicate branching. For example, the first term in (5) represents the event that a particle outside of an n-segment of 
odd parity jumps in, thus creating an n-segment of even parity. 

For the case of immediate reactions, r = 1, the processes indicated in (5) occur regardless of the state of target 
sites. It is easy to show (see [26, 27, 28, 29]) that: 

F n = Pr(r^^.) = ^[(1 - Gx) + (G n - G n+1 )} , (6a) 

H n = I'n^—. = i[(l - Gi) - (G n - G n+1 )} . (6b) 



Using these relations, Eq. (5) becomes 

d 



G n (t) = 2(r + n)(G B _i - 2G n + G n+1 ) . (7) 
at 
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r and f2 are the rates of hopping and branching, respectively, and the factor of 2 accounts for the events in (5) taking 
place also at the left edge of the interval. The case of n = 1 requires a special equation, since Go is undefined. Taking 
into account all the ways G\ might change, we find 

= 2r(l - 2G X + G 2 ) + 2n(G 2 - Gi) . (8) 
Thus, Eq. (7) may be understood to be valid also for n = 1, provided that one uses the boundary condition 

Additionally, since the G„ are probabilities, we have 

< G n (t) < 1 . (10) 

Eq. (7), with the boundary conditions (9), (10) may be analyzed exactly for a variety of initial conditions [27]. Here we 
merely observe that the only steady state solution supported by these equations is G n = l,n = l,2,..., corresponding 
to the absorbing state p s = 0. Indeed, the 2-BAW model with r = 1 lacks an active phase and exhibits no transition. 

IV. PARITY INTERVALS APPROXIMATION FOR THE GENERAL CASE 

Consider now the case of finite reaction probability, r < 1. A transition about some critical probability value r c is 
known to take place, from the absorbing state into an active (non-empty) phase [25], and our goal is to capture this 
transition, if only in an approximate fashion. We shall assume that r = O = 1, without loss of generality [35]. 

Consider the first process on the r.h.s. of (5). For r < 1, there is a difference in the reaction rate depending on 
whether the target site is empty or occupied: 

n < n— 1 < 7i—l i 

^^^m» = ^^^h • + r i ! •• . 

Unfortunately, events such as ^^^m o • and i i • • cannot be expressed in closed form in terms of the G n ■ We 

therefore rewrite the process in a way that the terms associated with r < 1, and which cannot be expressed in closed 
form, appear as a perturbation, proportional to s = 1 — r: 

n < n — 1 < n— 1 < 

n-1 , . n-1 . . 

= ^^^mo«+ i [ ••— (1 — r) I m (Ha) 

n n—1 

Likewise, the remainder of the diffusion events in (5) may be rewritten as 

n— 1 > n — 1 n—1 \ 

I I « = I n-.sr | » . (lib) 

(lie) 

(lid) 



n-1 



□ • — L 


!• 


— .S 1 


!•• 




n 


n-1 




□ • — C 




s ^^^m 


■ •• 




n-1 


n-1 




■ • = 1 




S 


■ •• 



The branching events too require close inspection of the target sites. For example, we rewrite the first branching 
event of (5) in a perturbative fashion: 

n >^t^ n-1 n-1 n-1 <sf^» n-1 

^^^^m» =^^o»o+ri itto|r^^MO«»|r i »»« 

, * n—1 n—1 n—1 n—1 n—1 , * 

= (1 — r) ^^^h • + r ^^^h • + r i m o + r ^^^m •• + r i m « (12a) 

n—1 - , n 

= S ^^^M • + (1 — -s) ^^^H* . 

The remainder of the branching terms are similarly expressed as 

n \ */\\t n 2 n 1 

i i » =s i ~ io»o + (1 — s) i — i « , (12b) 



n "ST* n-1 s n „ , 

i i « = s i i Qto + (1 — s) i i « , (12c) 

„_1 i/K „_2 ,, . n-1 

^^H* - s^HO.cMl (12(1) 
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So far everything is exact. In order to proceed, we approximate the terms proportional to s in the simplest possible 
way, by neglecting correlations. Thus, for example, we write 

Pr(r^=^..) « Pr(r^=^.) Pr(.) = F n ^(l - Gi) , (13) 

for the problematic term in (11a), and similar expressions for the ones in (llb)-(lld). For the problematic terms 
proportional to s in (12), we introduce the notation 




and 

Pr(»») = x , 

and approximate in the same spirit as above, 

Pr(^^o . o) w Pr(^^o) Pr(.o) = h n - lV , (14) 

and similarly for the other terms. Since y + z = G\, x + z = G 2 , and x + 2y + z = 1, it follows that x = §(1 — 2Gi +G2), 
y = ^(1 — G 2 ), and z = ^(— 1 + 2Gi + G2) are expressible in terms of the G„ in closed form. Since /„ + F n = G n and 
h n + H n = 1 — Gi, it follows that /„ = ^(— 1 + G\ + G n + G n+ i) and h n = 5(1 + Gi — G„ — G„ + i) are too given in 
terms of the G„ in closed form. 

We are now ready to write down a closed evolution equation for G n . Starting from Eq. (5), using the representations 
of Eqs. (11) and (12) with their respective approximations, collecting terms and rearranging, we obtain 

jG n = X -s{\ - G 2 )G„_ 2 + (2 - 3 S + 2sG 1 )G„_ 1 + i(-8 + 7s - AsG 1 + sG 2 )G n + (2 - s)G n+1 , (15) 

valid for n > 3. As expected, this reduces to (7), with T = Q = 1, in the limit s — > 0. The equation for n = 2 is 
exactly the same as (15), provided that one adopts the boundary condition 

Go = 1 . (16) 

Gi requires a separate evolution equation that we find by considering all the events that contribute to dG\/dt, and 
write down, in the spirit of Eq. (15), as 

d <- -» <- As, A\ 

— Gi = 2(»» + • — o» + •• — o» ) 

w 2{(1 - a)x + [(1 - GO - ra ] - y + (l - s ).t - [(1 - s )y + syd]} (17) 
= 2 - 6G1 + 4G 2 + s(-2 + 5d - 4G 2 + GiG 2 ) . 

To compute the steady state, we set dG n /dt = 0. Eq. (15) yields a recursion relation for the G„, with constant 
coefficients (that depend partly on Gi and G 2 ). Applying the ansatz G„ = A™ we find a cubic equation for A, with 
roots 

4 - 5s + AsGi - sG 2 ± yftZ -5s + 4sGi - sG 2 ) 2 + Ss(2 - a)(l - G 2 ) 
Ao = l, A ± = . (18) 

Thus G„ has the general solution 

G n = A\\ +B\ n _+C , 

where A, B, C are constants to be determined from boundary conditions. Since < s,Gi,G 2 < 1 (in the active 
phase) it follows that |A±| < 1, and G = lim^oo G n = Goo- Suppose that the initial distribution of particles is 
random, at density po, then G n (t = 0) = \ + |(1 — 2po)™, and Goo(0) = 1/2 [27]. Because our model conserves parity, 
it follows that C = 1/2. Furthermore, the boundary condition (16) implies B = | — A. The remaining coefficient, A, 
could be found from the relations 

d = AX+ + (1 - A)A_ + 1/2 , (19a) 



G 2 = A\\ + (1 - A)A 2 _ + 1/2 , (19b) 
G 2 - 2S - 2+ J 6 -t S)Gl ■ (19c) 
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FIG. 2: Order parameter p as a function of r as obtained from (a) Monte Carlo simulations (diamonds), (b) parity interval 
approximation (lowermost curve), and (c) one-site cluster approximation (uppermost curve): the curve shows the stable root 
(see figure 1), (observe that p corresponds to c s in the cluster approximation section). Note the absence of a transition point 
in the third case. 

The first two relations are required for self-consistence, since A± depend on G\ and G2 (as well as on s). Eq. (19c) is 
derived from (17), with dGi/dt = 0. Because Eqs. (19) are unwielding, in practice we set a numerical value for G\ 
and search for (A, s, G2) that satisfy (19), using Mathematica. We thus find a kinetic phase transition which we next 
compare to simulations. 

The critical point r c = 1 — s c , and the behavior of the order parameter in its vicinity, can be obtained analytically, 
by perturbing (19) about the values A = 0, G x = 1. We thus find r c = ± V33 - \ = 0.457427 ... and p s ~ (r c - r) 13 , 
with P = 1 . The value of (3 = 1 is in accordance with the expected from a mean-field like solution. 

We can also obtain analytically the behavior of the order parameter in the neighborhood of r = 0. The answer, 
c s ~ 1 — y/2r, is the same as the one obtained from the one-site and two-site cluster approximations, and seems to fit 
the numerical data from simulations exactly. Clearly, the role of fluctuations diminishes as r decreases away from the 
critical point, and it is conceivable that their impact is negligible in the limit of zero reaction. 

Comparison to Simulations 

In order to test the limit of validity of the above calculations, we have performed a Monte Carlo simulation of 
the 2-BAW model. We use a system size L = 10 4 , and implement the previously discussed rules. The initial 
configuration consists of a randomly half-occupied lattice; we have verified that the obtained long-time asymptotic 
results are insensitive to variations of the initial condition. Plotting the average density as a function of time in 
a double logarithmic plot, we identify the critical point with the value of r for which a separatrix, between curves 
that tend asymptotically to a constant (active phase) and curves that bend downwards (absorbing phase) converging 
asymptotically to a t -1 / 2 decay [36] is obtained. We thus find r c = 0.495(10). From the slope of the plot the density 
decay critical exponent 6 = 0.28(1) (extending for about four decades), is determined. It is in excellent agreement 
with previous measurements in the DP2 class [4]. The stationary density values are represented in Fig. 2. 

A point worth emphasizing is that finite-size corrections differ from those usually encountered in systems with 
absorbing states: Rather than finding a larger stationary order parameter for smaller systems, as is common, here 
one observes a faster decay to the absorbing phase. Owing to this the exponent ratio (3 fv cannot be determined from 
finite-size scaling analysis. 

We have also measured the order-parameter critical exponent, by plotting the average stationary density as a 
function of the distance to the critical point (Fig. 3). Our finding of j3 = 0.92(5) is again in good agreement with the 
commonly accepted value [4] . Other exponents can certainly be measured, but with 9 and {3 we can already guarantee 
that the model is indeed in the DP2 class. 

A more stringent test of the mean-field parity interval approximation is obtained by comparing its predictions for 
G n to the numerical results of Monte Carlo simulations, for different values of the control parameter r (Fig. 4). In all 
cases, we expect G n — > 1/2 for large n, since G-o = 1/2 and parity is conserved. The relevant question is how well the 
approximation captures the G n for n small. Because the typical distance between particles grows as one approaches 
criticality, we plot the results against v = pn rather than n, in order to compare them better. 
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FIG. 3: Log-log plot of the order parameter vs. the distance to the critical point. Shown are simulation results (curve with 
diamonds) compared to the prediction from the parity interval approximation (solid curve). Slopes of 1 (top, dashed line) and 
0.92 (bottom, dotted line) are shown for comparison. The latter represents the accepted value of the order parameter exponent 
/3 for the DP2 universality class in d = 1. Nearby criticality, the mean- field solution converges to (3 = 1. Surprisingly, the 
mean-field slope agrees quite well with the expected value of 0.92, when r c — r > 0.1. 




FIG. 4: G„ in the active phase for three different stationary activity values: (from the bottom to the top) 0.633, 0.2, and 0.1, as 
computed in a) Monte Carlo simulation for r = 0.082, r = 0.38 and r = 0.45 respectively (curves with symbols) and b) interval 
approximation (continuous lines). Instead of plotting as a function on n in order to make more meaningful the comparison we 
use v, defined as n divided the average interparticle distance in each case. 

For small values of r (as shown in Fig. 4), deep into the active phase, spatial correlations play a minimal role and 
the theoretical prediction matches simulations remarkably well, for all v. As r is raised toward the critical point, 
correlations play a larger role and the fit worsens. For example for r — 0.38 (p = 0.2) the theoretical curve fits 
experiments well only up to v m 2 (and, of course, also at v ^S> 1). Closer still to the critical point, at r = 0.45 
(p = 0.1), the theoretical prediction is good only up to v w 1. The large deviations observed for v > 1 in this last case 
indicate that the state of the system near the DP2 transition, despite being dilute, is highly ordered and correlated. 

The approximate curves, give converge to the exact value for G n — > 1/2 for asymptotically large values of v. 

V. DISCUSSION AND FUTURE PERSPECTIVES 

We have presented an approximation, based on the method of parity intervals, for the analysis of the DP2 transition 
observed in the 2-BAW model in one dimension with finite reaction rates. The key ingredient of our approach is 
that it respects the parity conservation implicit in the dynamical rules of the 2-BAW model: parity conservation is 
here responsible for the Zi symmetry that underlies transitions in the DP2 class, and provides a surprisingly good 
description of the transition at a mean-field level. 
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An interesting finding is the fact that spatial correlations play no role in the limit of zero reaction rate, r — > 0. In 
this "saturation" limit, the dependence of the order parameter upon the critical field, p ~ 1 — V2r, is captured even 
by the simplest one-site cluster approximation. 

Extensive numerical simulations of the DP2 transition suggest a value of [3 — 0.92(2) for the order parameter 
critical exponent [4]. On the other hand, an exact value (3 — 1 was conjectured by Jensen and later supported by 
numerical simulations combined with Pade approximants [16]. Regrettably, we have little to add to this controversy. 
Our mean-field results away from criticality (r c — r > 0.1) are consistent with (3 = 0.92, while f3 — > 1 as we approach 
the transition (Fig. 3). However, the large spatial correlations observed near criticality arc not faithfully modeled by 
our mean-field like approach (Fig. 4), and the values we obtain for r c — r < 0.1 are unreliable. It would be desirable 
to improve our approximation so as to better model regions closer to r c . Our attempts in this direction have been 
futile, so far. Arguably, the simplest improvement would be to replace the approximation of Eq. (13) by the more 
accurate: 

p , n-i . PrH^.)Pr(..) 
Pr(r=..)« ^ , 

and likewise elsewhere. However, on employing this scheme, instead of an improved result we find that the transition 
disappears. A similar phenomenon is known to occur also in mean-field cluster approximations, where sometimes 
increasing the cluster size yields lower quality predictions. 

Other interesting open prospects include using the parity interval approximation for the analysis of dynamical 
aspects of the DP2 transition. We have here systematically assumed that all time derivatives are zero, thereby 
accessing the steady state alone. Analyzing the very same equations with the full time dependence built in should 
yield a prediction for the order parameter time decay. Additionally, the equations could be modified to describe a 
spatially inhomogeneous system. In that case, one could study the dynamics of spreading. 
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